Cardiac mechanics Benchmark - Problem 2#
This code implements problem 2 in the Cardiac Mechanic Benchmark paper
Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.
from pathlib import Path
import dolfin
try:
from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
from dolfin import DirichletBC, Constant, interpolate, Mesh
import pulse
from fenics_plotly import plot
here = Path(__file__).absolute().parent
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 here = Path(__file__).absolute().parent
NameError: name '__file__' is not defined
geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["benchmark"])
# geometry = pulse.geometries.benchmark_ellipsoid_geometry()
2024-01-29 16:19:22,465 [925] INFO pulse.geometry_utils:
Load mesh from h5
# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["CC"] = 10.0
material_parameters["bf"] = 1.0
material_parameters["bfs"] = 1.0
material_parameters["bt"] = 1.0
material = pulse.Guccione(parameters=material_parameters)
# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
return DirichletBC(
V,
Constant((0.0, 0.0, 0.0)),
geometry.ffun,
geometry.markers["BASE"][0],
)
# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
# Solve problem
pulse.iterate.iterate(problem, lvp, 10.0, initial_number_of_steps=200)
2024-01-29 16:19:24,404 [925] INFO pulse.iterate: Iterating to target control (f_28)...
2024-01-29 16:19:24,404 [925] INFO pulse.iterate: Current control: f_28 = 0.000
2024-01-29 16:19:24,405 [925] INFO pulse.iterate: Target: 10.000
2024-01-29 16:19:24,555 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:24,556 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:24,556 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:24,557 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:19:25,239 [925] DEBUG UFL_LEGACY: Blocks of each mode:
99 full
2024-01-29 16:19:25,449 [925] DEBUG UFL_LEGACY: Blocks of each mode:
18 full
2024-01-29 16:20:25,159 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,160 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,160 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,161 [925] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-01-29 16:20:25,319 [925] DEBUG UFL_LEGACY: Blocks of each mode:
10 full
2024-01-29 16:20:25,515 [925] DEBUG UFL_LEGACY: Blocks of each mode:
3 full
*** Warning: PETSc SNES solver diverged in 2 iterations with divergence reason DIVERGED_FNORM_NAN.
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 44),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 82),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 107),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 140),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 173),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 255),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 296),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 337),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 378),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 419),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 460),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 501),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 542),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 583),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 608)],
[Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 42),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 80),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 105),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 138),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 171),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 253),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 294),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 335),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 376),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 417),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 458),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 499),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 540),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 581),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 606)])
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
endo_apex_marker = geometry.markers["APEX_ENDO"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)
print(
("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
endo_apex_pos[0],
),
)
Get longitudinal position of endocardial apex: 30.7991 mm
epi_apex_marker = geometry.markers["APEX_EPI"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)
print(
("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
epi_apex_pos[0],
),
)
Get longitudinal position of epicardial apex: 32.0417 mm
V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()